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Volume-weighted statistics of large-scale peculiar velocity is preferred by peculiar velocity cos¬ 
mology, since it is free of the uncertainties of galaxy density bias entangled in observed number 
density-weighted statistics. However, measuring the volume-weighted velocity statistics from galaxy 
(halo/simulation particle) velocity data is challenging. Therefore, the exploration of velocity assign¬ 
ment methods with well-controlled sampling artifacts is of great importance. For the first time, we 
apply the Kriging interpolation to obtain the volume-weighted velocity field. Kriging is a minimum 
variance estimator. It predicts the most likely velocity for each place based on the velocity at other 
places. We test the performance of Kriging quantified by the E-mode velocity power spectrum from 
simulations. Dependences on the variogram prior used in Kriging, the number nt of the nearby 
particles to interpolate, and the density np of the observed sample are investigated. First, we find 
that Kriging induces 1% and 3% systematics at fc ~ O.l/iMpc”^ when np ~ 6 x 10~^(/i~^Mpc)~^ 
and np ~ 6 x 10“^(h“^Mpc)“®, respectively. The deviation increases for decreasing np and in¬ 
creasing k. When np < 6 x 10“'^(h“^Mpc)“^, a smoothing effect dominates small scales, causing 
significant underestimation of the velocity power spectrum. Second, increasing Uk helps to recover 
small-scale power. However, for np < 6 x 10“"^(/i“^Mpc)“® cases, the recovery is limited. Finally, 
Kriging is more sensitive to the variogram prior for a lower sample density. The most straight¬ 
forward application of Kriging on the cosmic velocity field does not show obvious advantages over 
the nearest-particle method [Y. Zheng, P. Zhang, Y. Jing, W. Lin, and J. Pan, Phys. Rev. D 88, 
103510 (2013)] and could not be directly applied to cosmology so far. However, whether potential 
improvements may be achieved by more delicate versions of Kriging is worth further investigation. 

PACS numbers: 98.80.-k, 95.36.-l-x, 98.80.Bp 


I. INTRODUCTION 

The cosmic velocity field carries cosmological informa¬ 
tion. Peculiar velocity probes the growth rate of the Uni¬ 
verse. This makes it powerful to separate gravity and 
dark energy models (e.g.. Refs. p]-H^). Obtaining the 
observed number density-weighted and mass-weighted 
velocity is straightforward from galaxy velocity data and 
simulations. However, statistics involving the observed 
number density-weighted velocity suffers from uncertain¬ 
ties in galaxy density bias. On the contrary, the volume- 
weighted velocity statistics is free of uncertainties from 
galaxy bias. Thus, it is of particular importance since it 
is desired for the purpose of cosmology. However, mod¬ 
elling of the volume-weighted velocity is nontrivial due 
to the samp ling ar tifacts in velocity assignment methods 
(e.g.. Refs, prfil^ ). These sampling artifacts arise from 
assigning velocity on grids from highly nonuniformly dis¬ 
tributed particles/galaxies, especially for the populations 
with a low number density. 

These sampling artifacts in observation bias the veloc¬ 
ity power spectrum measured from galaxy velocity data. 
The same artifacts in simulation hamper the theoretical 
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understanding of the velocity field. A biased theoreti¬ 
cal understanding can lead to biased cosmological con¬ 
straints, even the velocity measurements are free from 
sampling artifacts like the one inferred from redshift- 
space distortion (RSD) measurement [1^. To proceed, 
two aspects could be improved. A straightforward one is 
choosing a reasonably good velocity assignment method 
with negligible sampling artifacts. Alternatively, under¬ 
standing and modelling the sampling artifacts for a cho¬ 
sen assignment method could correct the systematics and 
improve the velocity statistics measurement. 

In the literature, several methods have been adopted to 
estimate the volume-weighted velocity field. The Voronoi 
tessellation (VT) method is a zeroth-order interpolation 
scheme. One first constructs Voronoi tessellation from 
a set of nodes (i.e., particles/haloes) which forms the 
Voronoi polyhedral. Each Voronoi polyhedral contains 
only one particle and shares a wall with another poly¬ 
hedral. Inside a Voronoi polyhedral, the velocity is ap¬ 
proximated as a constant, which is just the velocity of 
the particle inside. Smoothing this space-filling velocity 
field, one could obtain the velocity on regular grids. The 
Delaunay tessellation (DT) method [2lj is a linear inter¬ 
polation scheme. One first constructs Delaunay tessella¬ 
tion, which is the dual of Voronoi tessellation and forms 
the Delaunay tetrahedron. The velocity gradient inside 
a Delaunay tetrahedron is approximated as a constant 
and is determined by the velocity of the four vertices. 



2 


Also, the volume-weighted quantities are determined by 
applying a smoothing on the interpolated velocity field. 
In Bernardeau and van de Weygaert’s work [22j |. the 
above two methods were proposed and found to agree 
well for a considerable range of situations. Bernardeau 
et al. [ 2 ^ found that the velocity divergence probabil¬ 
ity distribution function (PDF) measured by the above 
two methods also successfully produces the dependence 
on the matter density parameter. Pueblas and Scocci- 
marro adopted the DT method to define the velocity 
field, showing that the power spectra of velocity diver¬ 
gence and vorticity measured in this way are unbiased 
and have better noise properties than for standard inter¬ 
polation methods that deal with the mass-weighted ve¬ 
locities. However, theoretical modelling by Zhang et al. 
[Tsi l argues that none of these methods are free of sam¬ 
pling artifacts, in particular at the stringent 1% level 
required by Stage IV dark energy projects. 

In the works by Zheng et al. [I^ and Koda et al. [^ . 
the nearest-particle (NP) method was proposed and ap¬ 
plied. It is just the first step of the VT method that 
assigns the velocity of a regular grid point as the veloc¬ 
ity of the nearest particle to it. In this sense, it is the 
most straightforward velocity assignment method. With¬ 
out applying a smoothing, it alleviates artificial suppres¬ 
sion of small-scale random motions. Nevertheless, the NP 
method suffers from sampling artifacts, first detected by 
Zheng et al. [l^. More robust and systematic detection 
of the sampling artifacts from the NP method through 
N-body simulation is reported by Zheng et al. [H. It 
causes ^ 12% underestimation of the velocity power spec¬ 
trum at fc = O.l/iMpc”^ for samples with mean number 
density ^ 6 x 10“"‘(/i“^Mpc)“^. An advantage of the 
NP method is that the sampling artifacts are relatively 
straightforward to understand. Theoretical modelling of 
this sampling artifact helps to correct the systematics 
to some extent [l^, [l^. However, we are still not able 
to match the stringent requirement (1% modelling ac¬ 
curacy) of Stage IV dark energy projects. Further im¬ 
provements on modelling the sampling artifacts and/or 
better velocity assignment methods are desired. This pa¬ 
per presents our tests of Kriging as a velocity assignment 
method. 

Kriging, originated in geostatistics, is a method of in¬ 
terpolation for which the interpolated values are mod¬ 
elled by a Gaussian process governed by prior covari¬ 
ances. It can be understood as a linear prediction or 
a form of Bayesian inference. The Kriging method starts 
with a set of values observed and assumes a prior distri¬ 
bution over the whole field. A new value can be predicted 
at any new spatial location, by combining the Gaussian 
prior with a Gaussian likelihood function for each of the 
observed values. The resulting posterior distribution is 
also Gaussian, with a mean and covariance that can be 
simply computed from the observed values, their vari¬ 
ance, and the kernel matrix derived from the prior. As¬ 
suming the Gaussian velocity field, it gives the most likely 
velocity at the interpolated position based on the con¬ 


figuration of points with known velocity. It is accurate 
when the interpolated position overlaps with one of the 
observed points. The weighting derived from the Kriging 
method does not depend on the observed values. 

It is worth noting that the Kriging method is also use¬ 
ful in other branches of astrophysics. For example, in 
the field of weak gravitational lensing, the presence of 
the point spread function (PSF) biases the galaxy shape 
measurement. To correct it, one needs to interpolate the 
PSF form at the galaxy position based on the star im¬ 
ages nearby. It has been found that the Kriging method 
generally performs better than other methods, including 
polynomial interpolation, radial basis functions, and De¬ 
launay triangulation, in the context of weak lensing . 
We are therefore strongly motivated to test the perfor¬ 
mance of Kriging on cosmic velocity fields. 

This paper is organized as follows. In Sec. im we in¬ 
troduce the Kriging method in general cases and in the 
application on the cosmic velocity field. The simulation 
we use is briefly described in Sec. uni The statistics with 
which we are most concerned, the E-mode velocity power 
spectrum obtained from the Kriging method, is also in¬ 
troduced in this section. We present the dependence of 
the Kriging method performance on various factors in 
Sec. HVl We conclude and discuss in Sec. |Vl 

II. KRIGING INTERPOLATION 

Under suitable assumptions on the prior, Kriging gives 
the best linear unbiased prediction. The method is 
widely used in the domain of spatial analysis and com¬ 
puter experiments. The technique is also known as the 
Kolmogorov-Wiener prediction [271 l28l|. 

A. General idea of Kriging 

We assume that we want to estimate the field at a 
given position x* from Uk observed data points at 
with a value of yi = ?/(xi). Kriging looks for the value of 
the field at this position as a weighted linear combination 
of the nearby values at known positions, 

= 

i 

The weighting W is estimated such that the estimator is 
unbiased, and it minimizes the error with respect to the 
data according to the mean square variation. Therefore, 
Kriging relies on the variogram of the data as a prior to 
interpolate. The variogram is the mean-square variation 
of field values as a function of the separation, 

7(r) = ^ ([2/(x -f r) - ?/(x)]2) , (2) 

in which the ensemble average is over all spatial positions. 
Usually one assumes that the variogram only depends on 
the distance under the assumption of isotropy. 
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The weighting IV is solved from the following Kriging 
system: 
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Here, 'jij is the variogram matrix between observed 
points, with "fij = = 7 (|xi — Xj|). It has a dimen¬ 

sion of rifc X rifc. Unity vector 1^ is set for the unbiased 
condition '^^Wi = 1. is the variogram vector be¬ 
tween the interpolated point and the observed points, /r 
is Laplacian multiplier. We refer the readers to the Ap¬ 
pendix for the detailed derivation and some useful prop¬ 
erties of the Kriging method. 

One could measure this variogram from the observed 
data set itself numerically. For general usage, one could 
adopt a suitable model and fit the parameters from the 
data set to obtain an experimental variogram. Gaus¬ 
sian, spherical, and exponential variogram models are 
frequently used in applications. 

Usually the variogram function approaches to zero 
when r —>■ 0. For some cases, field values could change 
significantly even for small separation. One example is 
the existence of a measurement error for an observed sam¬ 
ple. In this case, one obtains different values for multiple 
observations even at the same point. Another example in 
cosmology is the nonlinearly evolved velocity field. After 
shell crossing, multistreams exist in cosmic flow, leading 
to a nonzero variogram at r = 0. This variogram behav¬ 
ior is called the nugget effect. With the nugget effect, 
Kriging always tries to smooth the field to suppress un¬ 
certainties existing within a small distance. Thus, in this 
case, Kriging interpolation no longer gives the observed 
value even that the position to interpolate approaches 
one of the observed points. 

B. For cosmic velocity field 

The cosmic velocity field at large scales (fc ~ 
O.lhMpc”^) is well approximated as curl free. Thus, the 
large-scale behavior is well described by the velocity di¬ 
vergence 9{x) = —V • v(x)/i7. Here, the normalization 
H is the Hubble parameter. The large-scale velocity is 
coherently evolved with the density field, which is ho¬ 
mogenous and isotropic. Thus, we adopt homogenous 
and isotropic velocity variogram in the Kriging interpo¬ 
lation. This variogram only depends on the separation 
of the two points: 

= i{nj) = i(|v(x,) - v(xj)p) . (4) 

Averaged over the directions of r^-, the correlation be¬ 
tween different spatial velocity components vanishes. 
Thus, each velocity component is Kriging interpolated 
from the same velocity component of other points. 

For general Kriging interpolation, a power-law fit var¬ 
iogram is sufficient to obtain a reasonable good estima¬ 
tion. However, a coarse fit of the velocity variogram 


may lead to an intolerant error in the estimated velocity 
power spectrum for precision cosmology. On the the¬ 
oretical aspect, we do not have an accurate variogram 
model for cosmic velocity due to the nonlinear evolution. 
However, we could obtain the linear velocity power spec¬ 
trum/correlation function prediction, and thus the linear 
variogram. Although it neglects the nonlinear effect that 
is difficult to model, the linear prediction is a good choice 
as the prior for the performance test. 

We do not include the nugget effect in this performance 
test for several reasons. First, although there exist ve¬ 
locity measurement errors for haloes/galaxies in observa¬ 
tion, we only deal with simulation particles with an ac¬ 
curate velocity measurement. Second, we only need the 
most likely velocity on grid points predicted by Kriging 
instead of the full posterior distribution. The full poste¬ 
rior distribution is useful in Kriging fitting, in which the 
Nugget is an important component. Finally, we want to 
keep the property that the velocity of the grid point is 
assigned to the velocity of the particle passing it. This 
will not suppress the random motions responsible for the 
fingers-of-god effect at a small scale where the particles 
are dense. 

Note that the velocity correlation (r) = 

(r!i(xi)t!j(x 2 )) between the ith velocity component 
at position xi and jth component at X2 = xi -|- r can be 
decomposed into two correlation functions tp± and 

(r) = V'-L {r)5ij + [V'li [r) - (r)] ^ . (5) 

Thus, the velocity correlation/variogram is anisotropic 
on all scales (also see Fig. 7 and the related descrip¬ 
tion by Zheng et al. 0). Considering this anisotropy in 
the velocity correlation/variogram, in principle, we could 
Kriging interpolate one velocity component from all the 
three velocity components of other points. This is an 
issue for future investigation. 

Toward small scales, the velocity variogram is inho¬ 
mogeneous for significant environment dependence. For 
example, the velocity difference is small along a filamen¬ 
tary structure and large in the perpendicular directions. 
This anisotropic variogram depends on the detailed en¬ 
vironment. Thus, for better reconstruction of the small- 
scale velocity field, we expect better performance with a 
position dependent anisotropic variogram. 

It is interesting to test the performance of Kriging with 
an inhomogeneous and anisotropic variogram. However, 
the main focus of this paper is the performance test of the 
most straightforward application of Kriging to measure 
the volume-weighted velocity on scale k ~ 0.1ft.Mpc“^ as 
the first step. We will discuss the possible extensions by 
adopting an inhomogeneous and anisotropic variogram in 
Sec. E and leave them for future work. 
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III. SIMULATION SPECIFICATION 
A. Simulation 

The simulation is run by a particlc-particle-particle- 
mesh code adopting standard flat Lambda cold dark mat¬ 
ter (ACDM) with rim = 0.268, Ha = 0.732, cts = 0.85, 
Us = I, and h = 0.71 (see Ref. 0). There are 
1024^ dark matter (DM) particles inside a box of size 
1200/i“^Mpc. We randomly sample the DM particles in 
simulation with fractions of / = 100%, 10%, 1%, 0.1% 
and 0.01% to see the dependence of the velocity statistics 
upon the sample density. The sample density np scales 
with / as 


For each case, we interpolate velocities on a 256^ uniform 
grid by the Kriging method from the nearest Uk = 200 
observed points. 

In observations, the volume-weighted velocity field is 
interpolated from halo/galaxy velocities. In principle, we 
should quantify the performance of Kriging on the halo 
field in simulation. However, due to the unknown halo 
velocity bias, we are lacking the correct volume-weighted 
halo velocity field even in simulation. While the num¬ 
ber density of the full DM sample in this simulation is 
sufficiently large and the velocity field constructed from 
the full sample could be treated as the right answer M- 
Thus, we choose to quantify the performance of Krig¬ 
ing on randomly selected DM particles. Although the 
randomly selected DM particles do not exactly share the 
same distribution as real haloes, this process keeps the 
clustering property of the matter in the Universe. 

Specifically, / = 0.1% and 1% correspond to num¬ 
ber densities of np ~ 6 x I0“^(h“^Mpc)“^ and 6 x 
10~^(/i~^Mpc)~^. A future wide spectroscopic survey 
like Euclid [3l| and DESl[32| could measure velocity data 
with this number density. Thus, the performance of Krig¬ 
ing for these two cases is of particular interest. 


the cosmic flow, producing the B-mode component on 
small scales [2^. Both components are related to the 
matter growth rate of the universe, thus carrying cosmo¬ 
logical information. 

Stage IV dark energy projects, such as MS-DESI (Big- 
BOSS), Euclid, SKA, and WEIRST, can measure the 
volume-weighted velocity power spectrum through RSD 
with 0(1%) statistical precision (e.g.. Ref. H). RSD is 
induced by a peculiar velocity. Thus, for RSD modelling, 
it is of crucial importance to understand the peculiar ve¬ 
locity field. An appropriate velocity decomposition has 
the potential to simplify and improve the RSD modelling 
[ 3 ^ . As the first step, the velocity is decomposed into an 
E-mode and B-mode component. The E-mode compo¬ 
nent is responsible for the large-scale Kaiser effect [l|. 

In this study, we focus on the E-mode component, 
which contains most of the cosmological information. 
First, the volume-weighted velocity field is Kriging in¬ 
terpolated on regular grids. Then, the interpolated ve¬ 
locity field is decomposed into an E-mode and B-mode 
component. 


V, = ve -I- vb , 


( 7 ) 


in which the E-mode component is curl free (V x ve = 0) 
and the B-mode component is divergence free (V • vb = 
0). This decomposition could be done using a fast Eourier 
Transform, 


VE 


(v* • k)k 


Vb = V* - Ve . 


( 8 ) 

(9) 


We measure the E-mode velocity power spectrum 


FVv(fc) = (veve*) , (10) 

and usually we present the result in the form of = 

fc3Rvv(fc)/27r2. 

We also obtain the same statistics for the velocity field 
obtained by the NP method. According to the conver¬ 
gence test by Zheng et al. [H , we consider the result 
from the NP method with / = 100% as a reference. 


B. Statistics 

Any vector field can be decomposed into an irrotational 
(gradient) part and a rotational (curl) part. Analogous 
to the electric and magnetic fields, we denote the first one 
with a subscript “E” and the later one with a subscript 
“B” [ 3 ^ . Decomposing the cosmic velocity field into the 
irrotational E-mode component and rotational B-inode 
component is widely used in cosmology studies. 

During linear evolution, the cosmic velocity field only 
has an E-mode component. Thus, at the early time or 
on the large scales, the velocity field is relatively easy 
to model. It could be completely described by its diver¬ 
gence, which is related to the matter density field [ 2 ^. 
Late-time nonlinear evolution induces shell crossing in 


IV. PERFORMANCE OF KRIGING 

The performance of Kriging depends on the variogram 
prior 7 (r), the number n^, of the observed points used 
to interpolate, and the density of the observed sample 
described by the sampling fraction /. Notice that the 
sampling fraction / (or equivalently, the sample density 
np) is not an input parameter in the Kriging method. 
When dealing with all dark matter particles in our sim¬ 
ulation, / = 100%. For massive and small haloes in ob¬ 
servation, the number density corresponds to / ~ 0.1% 
and ^ 1%, respectively. It is a fixed number when the 
observed sample is determined. However, it has a great 
impact on the performance of Kriging. Thus, we present 
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FIG. 1. Theoretically predicted variogram for the cosmolog¬ 
ical parameter used in simulation (solid line) and the vari¬ 
ogram from a Qm = 0.3 and i}rn = 0.236 flat universe (dotted 
line and dashed line, respectively). 


the dependence on / thorough the whole performance 
test. 

First, we quantify the performance of Kriging in the 
case that a theoretically predicted variogram is taken as 
the prior. Then, the dependence on rik is shown as a con¬ 
vergence test. Finally, the sensitivity on the variogram 
prior is tested by adopting inconsistent variograms. 


A. Theoretically predicted variogram 


We use the CLASS code [S^ to produce the velocity di¬ 
vergence power spectrum Agg(fc). The correlation func¬ 
tion is the Fourier pair of the power spectrum. 


,w = / 


dlnfc- 


A^g(/c) sin(fcr) 




kr 


( 11 ) 


Since the absolute amplitude of the variogram does not 
influence Kriging interpolation (see Appendix. lA 3() . we 
set the variogram as 


7 (r) = 1 - Cvv(?’)/fvv(0) . (12) 


This relation holds only if (v(a;)) = 0, true for the veloc¬ 
ity. The theoretically predicted variogram is presented 
in Fig. [T]as a solid line. It has a slope of zero beyond the 
velocity correlation length r ^ 150ft.“^Mpc. The most in¬ 
ner part of the variogram has slope of 2 by the definition 
of the correlation function, 

Cvv(?') = J dln/cPvv(fc)^^^^^^ 

= J d\nkP^^{k){l - ^{kr)'^ + 0{{kr)'^)) (13) 
= fvv(0)-y y dlnfcPvv(fc)fc^ H -, 

in which Taylor expansion is adopted in the second line. 


The right panel of Fig. [J is just the same result as 
Fig. 1 by Zheng et al. [l^ . For the NP method, it 
causes a ~ 12 % underestimation of the velocity power 
spectrum at k = O.l/iMpc”^ for the / = 0.1% case. This 
systematic underestimation increases with decreasing / 
and increasing k. However, also reported is an anomalous 
behavior, which is not expected from simple modelling 
of the sampling artifacts for the NP method. For the 
/ = 0 . 01 % case, the decreasing power turns up at small 
scales (fc > 0.08/iMpc“^). Two possibilities causing this 
anomaly are discussed in the appendix of Zheng et al. 
0- We refer the readers to Zhang et al. 0 and Zheng 
et al. [l^ for the detailed sampling artifacts study and 
modelling for the NP method. 

Adopting the theoretically predicted variogram, the E- 
mode power spectrum obtained from the Kriging method 
with Uk = 200 is presented in the left panel of Fig. [2] 
We find a well-reconstructed power spectrum for / = 
100%, 10%, and 1%. Obvious deviation only appears 
for very low sampling fractions / = 0 . 1 % and 0 . 01 %. 
For / = 0.1%, the power spectrum is suppressed at A: > 
O.l/iMpc”^, while for / = 0.01%, the suppression begins 
at a much larger scale (fc > 0.02/iMpc“^). The success for 
/ = 100 % to 1 % implies that this theoretically predicted 
variogram is reasonably close to the true one, although 
it neglects the nonlinear evolution of the velocity. The 
sudden failure for / = 0 . 1 % and 0 . 01 % shows that the 
sampling fraction has a great impact on the performance 
of Kriging. 

The mean particle separation Lp scales with / as 

Lp = /j-iMpc . (14) 

1024-' ^ ^ ' 

For / = 0.1% and 0.01%, the mean particle separation is 
11.72ft,“^Mpc and 25.25/i“^Mpc, respectively. The mean 
distance of one grid point to its nearest particle mea¬ 
sured in simulation is 7.77ft.“^Mpc and 15.55ft.“^Mpc, 
respectively. From Fig. [U we can see that for these 
separations of the / = 0 . 1 % case the variogram has val¬ 
ues of > 0.35 and 7 ,* > 0.2. This means that the 
velocity difference among nearby particles and between 
the grid point and nearby particles is large. Thus, in 
these situations, all the nearby particles obtain compa¬ 
rable weightings from the Kriging system, leading to a 
large smoothing effect. 

We define a characteristic length as 


L 


V 


A(V-v)^) V" 
1 (v^) ) 


(15) 


This characteristic length indicates on which scale the ve¬ 
locity varies significantly. If the mean particle separation 
Lp is larger than this characteristic scale Ly, the veloc¬ 
ity variation is not well sampled. Ly measured from the 
velocity field obtained by the NP method with / = 100% 
is Ly ^ 6.64h“^Mpc. For / = 100%, 10% and 1%, 
Lp = 1.17, 2.52 and 5.44/i“^Mpc, smaller than Ly. For 
these cases, the observed sample is sufficiently dense, and 
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FIG. 2. E-mode velocity power spectrum from the Kriging method is presented in the left panel in which Uk = 200 is adopted. 
The variogram prior in the Kriging method is predicted by the CLASS code. The one from the NP method is presented in the 
right panel. 




k (h/Mpc) k (h/Mpc) 


FIG. 3. Ratio of E-mode velocity power spectrum from the Kriging method with = 200 (left panel) and the NP method 
(right panel) to the reference one. The vertical dotted lines indicate the most concerned scale k = 0.1hMpc“^. Notice that for 
/ = 100% Kriging and NP agree with each other well within 1% at fc < O.l/iMpc”^ (left panel). This confirms the previous 
conclusion that the NP method is robust at fc < O.lfcMpc”^ when np > 0.6(/i~^Mpc)“®. 


Kriging has good performance, For / = 0.1% and 0.01%, 
Lp > Ly. In these situations, the velocity field is not well 
sampled, and therefore a smoothing effect dominates. 

We present the ratio of the E-mode velocity power 
spectrum from the Kriging method and NP method to 
the reference one in the left and right panels of Fig. 
According to the modelling of the sampling artifacts for 
the NP method by Zhang et al. [l^, the reference case 
only suffers from negligible sampling artifacts at scales 
k < 0.2/iMpc“^ (see Fig.3 by Zhang et al. [3). The sim¬ 
ulation results in Zheng et al. show that the sampling 
artifacts of / = 10% are only 1% at fc = O.l/iMpc”^. This 
also implies negligible sampling artifacts for / = 100%. 
Thus, the reference case is well approximated as the cor¬ 
rect answer at the scales we consider. The right panel of 
Fig. [3] is the same result as Fig. 2 by Zheng et al. [T^. 

We find that at scales larger than the smoothing dom¬ 
inating region power from the Kriging method is slightly 


larger than the reference one. The nonlinear velocity 
power spectrum is predicted and found to be lower than 
the linear prediction at small scal es (/K-i R ef. [3^ for a 
perturbative prediction and Refs. [24l.l38ll39l| for the sim¬ 
ulation result). Thus, we expect the true variogram to 
be lower than the linear prediction at small scales. The 
deviation of the theoretically predicted variogram from 
the true one may cause this overestimation in the power 
spectrum of the Kriging interpolated velocity field. 

Compared to the NP method, the most straightfor¬ 
ward Kriging does not show obvious advantages. For 
/ > 1%, the NP method underestimates the small-scale 
power more with decreasing /. This is an expected be¬ 
havior since for lower / the NP method assigns the same 
velocity to more grids, leading to a larger smoothing ef¬ 
fect, while for the Kriging method, the measured power 
is overestimated at intermediate scales and suppressed 
at small scales. Considering that the overestimation de- 
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pends on the variogram prior and we only use a linear 
prediction, a detailed comparison is less meaningful. In 
fact, in the following analysis (Sec. IIV C() . we find that 
the power amplitude depends on the cosmological param¬ 
eters used to predict linear variogram prior. This implies 
that the inaccuracy of the variogram prior could affect 
the performance of Kriging at the percent level. 

For / = 0.1% and 0.01%, a heavy smoothing effect 
dominates the performance of Kriging, while for the NP 
method, these two situations also produce anomalous be¬ 
havior, which is difficult to model. We emphasize that 
the sudden failure in Kriging and anomalous behavior 
in the NP method is physical. For such low sampling 
faction, the velocity field is not well sampled, and we do 
not expect good performance for any velocity assignment 
method. 

Notice that for / > 1% the sampling artifacts of the 
NP method are relatively straightforward to model and 
to some extent it could be corrected. Although the sam¬ 
pling artifacts are difficult to model for the more com¬ 
plicated Kriging method, we expect improvement by ap¬ 
plying a more delicate version of Kriging. We discuss the 
potential improvement that could be made in Sec. IVland 
leave this for a futrue investigation. 


B. Dependence on Uk 

In Fig. [U we can see that the variogram approaches to 
unity at ~ 150/i“^Mpc. This implies that particles with 
distance larger than ^ 150/i“^Mpc do not contribute to 
the prediction of the velocity on a given grid point. Thus, 
in principle, for each grid point, we only need to include 
all the particles with a distance less than ^ 150h“^Mpc. 
However, for high sampling fraction cases, this implemen¬ 
tation is still computationally expensive. Closer particles 
generally have a larger contribution to the estimation of 
the velocity on the grid point. Thus, we naturally choose 
to use the nearest rik particles from the grid point to 
Kriging interpolate. Which nu is suitable relies on many 
factors. 

One can image an extreme case with infinite particles 
in simulation. For any grid point, there always exists a 
particle passing it. Since Kriging is accurate in this case, 
rifc = 1 (i.e., reduced to the NP method) is enough. For 
general cases, whether the outer-part particles are im¬ 
portant depends on the steepness of variogram. A very 
steep inner part of the variogram means that the veloc¬ 
ity on the grid is very likely to be similar as close par¬ 
ticles. Thus, the weightings for outer-part particles are 
very small. In this case, excluding a part of the outer 
particles will not influence the interpolation result. 

The performance test with increasing Uk will greatly 
help us understand the Kriging method. Figure 0] 
presents the dependence on the particle number Uk used 
to Kriging interpolate. In the left plots, increases 
from 1 to 8, and in the right plots, it increases from 10 to 
160. The ratio of the interpolated velocity power to the 


reference one is presented to better show the dependence. 

For / = 100%, all the lines overlap with each other 
when k < O.lhMpc”^. This implies that for such a 
high density of particles the velocity on the grid could 
be determined well by only a few nearby particles. This 
result is consistent with the conclusion by Zheng et al. 

that for / = 100% the NP method is accurate. For 
/ = 10% and 1%, underestimation at small scales could 
be seen when only a few particles are used in interpola¬ 
tion. However, further increasing Uk increases the small- 
scale power toward (even over) the reference one. This 
behavior could be understood as follows. Inside a low- 
density region, several grid points share the same near¬ 
est particle, obtaining the same velocity from the NP 
method. Using only a few particles to Kriging interpolate 
also suffers from a similar sampling artifact. This sam¬ 
pling artifact assigns a similar velocity to several nearby 
grid points, therefore suppressing the small-scale power. 
Including more particles to interpolate gathers more in¬ 
formation. This helps to recover the difference between 
the nearby grid points, leading to the recovery of small- 
scale power. 

For / = 1%, 0.1%, and 0.01%, increasing Uk from 1 
to 8 gradually suppresses the power more. This reflects 
the smoothing effect of the Kriging method. However, 
further increasing from 10 to 160 recovers the small- 
scale power for / = 1% and increases the power at 0.02 < 
k < 0.2/iMpc“^ for / = 0.1% and 0.01 < k < O.lhMpc”^ 
for / = 0.01%. This also reflects the help by including 
more particles to interpolate. 

For / = 10% and 1%, increasing Uk to 160 almost 
fully recovers the power at all scales we consider. But 
for / = 0.1%, the power at scales k > O.l/iMpc”^ is still 
underestimated even rifc = 200. Further increasing tik 
will not help since new added outer-part particles are far 
from the grid point. Thus, this small-scale power is lost in 
the Kriging method. For / = 0.01%, the situation is even 
worse. The power is underestimated at fc > O.OlShMpc”^ 
when all particles are used. 

For / < 1% cases, increasing nu from 80 to 160 does 
not change the power significantly. Thus, the perfor¬ 
mance of Kriging converges at Uk ~ 100. 

Notice that the observed velocity sample has spatial 
clustering. Particles with coherent motion inside a small 
dense region do not provide independent information to 
Kriging system. Although this will not influence the re¬ 
sult, it reduces the performance of Kriging in the aspect 
of a waste of computational resources. For low sampling 
fraction cases like the halo population in our Universe, 
the mean separation is large, and we do not expect a large 
degradation. For a high sampling fraction, a preselected 
sample with fairer spatial distribution could improve the 
efficiency of Kriging interpolation. 

For the higher sampling fraction case, smaller rik is 
needed to reconstruct the power spectrum well. Thus, we 
could expect that for denser regions smaller Uk is needed. 
Determining Uk for each grid point from its density could 
also improve the computational efficiency. 





k (h/Mpc) 


FIG. 4. Ratio of the E-mode velocity power spectrum from the Kriging method to the reference one. The number of particles 
to Kriging interpolate increases from n* = 1 to 8 (left panel) and Uk = 10 to 160 (right panel). From top to bottom, the 
sampling fraction / decreases from 100% to 0.01%. The vertical dotted lines indicate the most concerned scale k — O.l/iMpc”^. 
Note that the y-axis scaling is different to better show the results. 


C. Sensitivity on the variogram prior 

For performance test on simulation, we know the true 
cosmology and could obtain a consistent theoretical var¬ 
iogram. When applying Kriging to observations, we do 
not know the cosmological parameters. Thus, we quan¬ 


tify the sensitivity on the variogram prior by adopting 
two inconsistent variograms in Kriging. These two var- 
iograms are produced by modifying 17^ from 0.268 to 
0.3 and 0.236 while keeping the universe flat, since the 
structure formation is very sensitive to ilm- 

Two inconsistent variograms we adopted are presented 
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FIG. 5. E-mode velocity power spectrum from the Kriging method using a consistent prior (solid line), and inconsistent ones 
(dotted and dashed line) are presented in the left panel. The ratio to the fiducial one is presented in the right panel. The 
vertical dotted line indicates the most concerned scale k = 0.1/iMpc“^. 
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in Fig. [H For larger (smaller) ilm, the (scaled) vari- 
ogram is slightly higher (lower) at the scales below the 
velocity correlation length. The E-mode velocity power 
spectrum obtained by the Kriging method with these in¬ 
consistent variograms is presented in the left panel of 
Fig. [SI and the ratio to the fiducial one is presented in 
the right panel. We find that the deviation amplitude is 
larger with decreasing /. An inconsistent variogram with 
larger (smaller) 17^ produces an underestimated (overes¬ 
timated) power spectrum. The underestimation (overes¬ 
timation) reaches 7—8% for / = 0.01% at fc ^ O.l/iMpc”^ 
and 5% for / = 0.1% at fc ~ O.S/iMpc”^. 

Notice that the power at the largest scale is insensi¬ 
tive to the variogram prior. It implies that Kriging keeps 
the consistency between the linear evolution of the den¬ 
sity and velocity field on large scales. For decreasing /, 
the deviation appears at larger scales with a larger am¬ 
plitude. Thus, for lower sampling fractions, Kriging is 
more sensitive to the variogram prior. 

The measured velocity power spectrum depends on the 
variogram prior, which itself depends on the cosmological 
parameters. This implies that the Kriging interpolated 
velocity field from a sparse halo population could not be 
used to constrain cosmological parameters. Thus, so far, 
we constrain the application of the Kriging method in 
the theoretical understanding of the peculiar velocity in 
simulation. 

With good understanding of the peculiar velocity field 
and a good enough velocity variogram/correlation func¬ 
tion model, one could try to fit the model parameter 
from the Kriging interpolated velocity field iteratively. 
This process could produce a volume-weighted velocity 
field with a self-consistent variogram, enabling the appli¬ 
cation for cosmology. We tried this iterative method by 
adopting coarse variogram models. We found the process 
is unstable since the resulting power depends significantly 
on the variogram model, especially on small scales, which 
is hard to accurately model. This hampers the applica¬ 


tion of Kriging in cosmology, unless we could obtain a 
sufficiently good variogram model. 


V. CONCLUSION AND DISCUSSION 

For the first time, we applied the Kriging method to 
obtain the volume-weighted velocity field from N-body 
simulations. We tested the performance of Kriging on 
the variogram prior 7 (r), the number Uk of the nearby 
particles to interpolate, and the density of the observed 
sample described by the sampling fraction /. 

First, we theoretically predicted a linear variogram 
consistent with the simulation. Adopting this theoret¬ 
ically predicted variogram, Kriging reconstructed the E- 
mode velocity power spectrum well for / > 1% at all the 
scales we consider. However, for / = 0.1% and 0.01%, 
a smoothing effect dominated small-scale power suppres¬ 
sion. This is a natural behavior when an observed sample 
is so sparse that the mean particle separation Lp is larger 
than the characteristic scale Second, we tested the 
rife dependence by increasing rife from 1 to 160. We found 
that increasing rife from 1 to 8 generally smoothed the 
small-scale power first, while further increasing rife from 
10 to 160 recovered the small-scale power to some extent. 
For / = 10% and 1%, the suppressed power could be re¬ 
covered at all the scales we considered, while for lower 
/, the power was eventually lost even when more par¬ 
ticles were included. The Kriging system converged for 
rife ^ 100. Finally, we adopted inconsistent variograms to 
check the sensitivity on the variogram prior. An incon¬ 
sistent variogram biased the velocity power spectrum at 
small scales. The deviation was larger for lower / and be¬ 
gan to appear at larger scales. Adopting the variogram 
predicted by a 0^ = 0.3 (0.236) flat universe, the un¬ 
derestimation (overestimation) of power reached 5% at 
k ~ O.lhMpc”^ for / = 0.1%. 

For / > 1% and / = 0.1% at k < 0.1hMpc“^, both the 
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most straightforward Kriging method and NP method 
generally performed at comparable levels. For lower / 
cases, both methods could not reconstruct the velocity 
field well. The Kriging method suddenly failed for these 
lower / cases, and the power suffered a heavy smoothing 
effect, while the NP method suffered from an unexpected 
behavior from simple sampling artifacts modelling, and 
was thus difficult to correct. 

There are several potential improvements to be made 
for using the Kriging method to obtain the volume- 
weighted cosmic velocity field in the future. First, as 
stated in Sec. IIVBI a better strategy for selecting the 
observed sample and determining the number of points 
used to interpolate could improve the computational ef¬ 
ficiency. Second, as discussed in Sec. IIIBl due to the 
fact that velocity field is a vector field, the correlation 
function/variogram is anisotropic on all scales. Adopt¬ 
ing the linear predicted but anisotropic variogram has 
the potential to improve the volume-weighted velocity 
reconstruction. In fact, after some calculations, we did 
find improvement by considering the anisotropy in the 
variogram. We expect this improvement to the Krig¬ 
ing method to be the most straightforward and potential 
one. Third, to better reconstruct the velocity on small 
scales {k > O.l/iMpc”^), taking the position-dependent 
anisotropic variogram into account is promising. To do 
this, we need understanding and modelling on the veloc¬ 
ity for different cosmic structures (void, sheet, filament, 
and cluster). Position-dependent Kriging interpolation 
is expected to suppress the smoothing effect for partially 
collapsed structures, leading to better small-scale recon¬ 
struction. Fourth, taking the correlation between all the 
velocity components into account, we could extend the 
Kriging method to interpolate each velocity component 
from all the velocity components. In fact, under the ap¬ 
proximation that the velocity is a gradient field, one can 
reconstruct one of the velocity components from another 
component equally as well as from the same component. 
Thus, whether this extension could better reconstruct the 
velocity field is worth exploration. Finally, we adopted 
the velocity variogram from the linear prediction, ne¬ 
glecting the nonlinear effect that is difficult to model. 
More knowledge on the velocity variogram (or equiva¬ 
lently, the velocity correlation function) would help us 
to better reconstruct the cosmic velocity field. Most im¬ 
portant of all, good modelling of the velocity variogram, 
especially on small scales, would enable us to fit a con¬ 
sistent variogram from the interpolated field iteratively. 
This would avoid the uncertainties induced by choosing 
the variogram prior. 

Better understanding of the cosmic velocity field will 
improve the Kriging interpolation. Conversely, to better 
understand the cosmic velocity field, one needs a good 
velocity assignment method with well-known or well- 
controlled sampling artifacts. For example, protohalo 
statistics predicts the physical halo velocity bias (&„ < I) 
(e.g.. Refs, However, a profound assumption 

in peculiar velocity cosmology is = 1 at a sufficiently 


large scale. Because of the severe sampling artifacts in 
measuring the volume-weighted velocity power spectrum 
for sparse populations, an accurate halo velocity bias 
measurement is hampered [13, IH, 113 • After appropriate 
corrections of the sampling artifact in the NP method, 
Zheng et al. Ea verified = I at fc < O.lhMpc ^ and 
2 = 0 — 2 for haloes of mass ^ 10^^ — 10 ^^Mq//i. Never¬ 
theless, the measured velocity bias is only accurate to a 
few percent due to the residual sampling artifact. Better 
correction of sampling artifacts is needed to measure 
with 1% accuracy. 

In this paper, we only quantified the performance of 
Kriging through the E-mode velocity power spectrum. 
The sampling artifact also affects other velocity statistics. 
For example, Zhang et al. [isj ] proved in the Appendix 
that the sampling artifact also affects the density-velocity 
correlation measured by the NP method. A similar sam¬ 
pling artifact is also expected in that measured by the 
DT method [13] ■ Although for the dense DM population 
that Jennings and Jennings [13| dealt with the sampling 
artifact is negligible, it must be appropriately corrected 
for a sparse halo population. Thus, the performance of 
Kriging on other statistics (e.g., the divergence, velocity 
shear, and B-mode component on small scales) is also in¬ 
teresting to explore. To study the B-mode velocity com¬ 
ponent, high-resolution analysis is essential since the B 
mode is significant on small scales and the correlation 
length is much shorter. To what extent Kriging could 
help us to better understand the cosmic velocity field is 
worth further investigation. 
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Appendix A: Kriging interpolation 


1. Derivation 

We assume that we want to estimate the field at a 
given position x* from rik observed data points at 
with value of yi = Kriging looks for the value of 

the field at this position as a weighted linear combination 
of the nearby values at known positions, 

y{y^*) = ^Wiy^ . (A.l) 

i 

The error in estimating y(x*) is e(x*) = y(x*) — ?/(x*). 
Since the estimator is modelled by a Gaussian process, 
to ensure that the model is unbiased, the unbiased con¬ 
straint is observed. 


£;(e(x,)) = 0 , (A.2) 

which leads to the unbiased condition • W = 1 . 

To find an estimator with minimum variance, we need 
to minimize £ 1 ( 6 ^(x*)): 

£i(e^(x*)) = • Cou(x,jXj) • W — Cou’^(xiX*) ■ W 

— W'^ ■ Cov{xiX^,) + Var{x^,) . 

(A.3) 


Under the assumption of isotropy, the elements in 
Cov{xiXj) are the correlation among the observed points 
f(|xi — Xjl), while the elements in C'ou(xiX*) are the 
correlation among observed points and the interpolated 
point .^(|xi - x*|). 

By the Lagrange multiplier method, solving this opti¬ 
mization problem results in the Kriging system 


w 


Cov{xiXj) 1 

-1 

Cov{'x.iX^,) 



1 ^ 0 


1 


Adopting the relation between the variogram and corre¬ 
lation 7 ( 7 -) = ^( 0 ) — ^{r), usually it is expressed in terms 
of the variogram: 


W 


lij 

1 

-1 

7« 




0 


1 


(A.5) 


Here the elements in 7 ^^- are the variogram among the 
observed points 7 (|xi — Xj|), while the elements in 7 ,* 
are the variogram among the observed points and the 
interpolated point 7 (|xi — x*|). 


In this case, it reduces to the NP method, i.e., setting the 
velocity on grid as the velocity of the nearest particle. 
Considering = 2, we have 


Wi 


0 

712 

I ■ 

-1 

71 * 

W2 

= 

712 

0 

I 


72 * 

^ . 


I 

I 

0 


I 




- 71 * + 712 
7l2 

- 72 * + 712 

712 


(A.8) 

(A.9) 


In the case of the grid point approaching one of the par¬ 
ticles (for example, particle 1 ), we have 71 * 711 = 0 , 

72 * ^ 712 , and [wi,W 2 \ —>■ [1,0]. Thus, if very close 
to one of these particles, the grid point will just obtain 
the velocity of that particle. In other words, the velocity 
assignment is accurate in this case, and the other parti¬ 
cles are not important. This is also true for any rife. It 
explains the result for the / = 100% case in which the 
power spectrum is insensitive to rife (see Sec. IIVBI) . 


3. Independence on the absolute amplitude of the 
variogram 

The velocity interpolated from the Kriging method 
only depends on the shape of the variogram (relative 
amplitude). Considering multiplying the variogram by 
a factor of a, the new weighting W and Laplacian mul¬ 
tiplier ft are solved from the new Kriging system 


a-Yzj 

1 


w 


Q!7i* 

1^ 

0 


A 


1 


The first line reads 

jijW + -fil = . (A.11) 

a 

The second line is the unbiased condition 

= (a.12) 

i 

From the above two equations one could easily find that 
the solution to this new Kriging system, Eq. (jA.lOl) . 
is the same as the solution to Eq. ©; i.e., W = W. 
The only difference is the irrelevant Laplacian multiplier 
fi = ay,. 


2. Trivial cases 

If only the nearest particle is used in Kriging, i.e., rife = 
I, we have 


w 


0 

I 

-1 

71 * 


1 



1 

0 


1 


* 


4. Singular case 

A general usage of Kriging interpolation is adopting 
a variogram fit from the same data used to interpolate. 
For applications without a high-precision requirement, a 
power law fitting 7 = ar^ is sufficient to obtain a rea¬ 
sonable good interpolated value. However, for a strong 
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linear case with 7 ^ = oc this Kriging 

system is degenerated. The proof is as follows: 


det(A) = 


1 


Id -rj|2 1 

1^ 0 


1 ^ 0 


r2 -hr 

M. t 1 J. 

^ - 2ri ■ Tj 1 


-2d • r? 1 


0 

h 

1—1 


1^ 0 


(A.13) 


Vi ■ rj 1 

-il^ 0 


= (- 2 )’"'= 

For rifc > 3, Ti = {xi,yi,Zi), we have 



xi yi zi 0 ^ 


Xi ■ 

Xuk 

_ 



yi ■ 

Vnk 

1 

r\'T' 


Zi • 

^rik 


Xuk Vuk 0 ^ 


0 • 

■ 0 


= 0 . 


(A.14) 


Thus, the right-bottom cofactor does not affect the de¬ 
terminant. We just set it as 0.5, 


det(A) oc 


1 -I- Ti • 0 


D-r^ 1 


l+Vi-Xj 1 

1- iT 1 

2- *- 2 


0^ i 

" 2 


= 1 -I- Ti ■ r,- 


1 xi yi zi 0^ 


1 Xuf, y-rik ^rik 0 


1 

Xi 
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-yTTk 


(A.15) 
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